%%%%%%%%%%%%%%
% fit_fj_gmm.m
%%%%%%%%%%%%%%

% Reset RNG
rand('state', 0);

% Data comes in as d x N matrix
% (i.e - rows are dimensions, columns are examples)

[d,N] = size(Data);

% Save obs data
eval(sprintf('save obs_%d_%d;', d, N));

switch lower(Covar_Type)
  case 'spherical'
    % disp('mixtures4.m does not handle spherical covariance matrices, defaulting to diagonal.');
    covoption = 1;
  case 'diagonal'
    covoption = 1; 
  case 'full'
    covoption = 0; 
  otherwise
    disp(sprintf('Unknown covariance type %s', Covar_Type));
    exit(1);
end

% set min & max no of components to try, and the threshold value
nmin = 1;
nmax = 10;
th = 1e-4;


% We need to run the F&J algorithm until it returns a sensible result
% Very occasionally it can return a GMM with no components, in which
% case re-running will generally fix it due to the random initialisation.

Num_Centres = 0;

% Use the IGDS speeded up version - mixtures4b
[Num_Centres, Priors, Means, Covariances] = mixtures4b(Data, nmin, nmax, 0, th, covoption);

if Num_Centres == 0
  cd /home/nph/Research/Matlab/Netlab; %matlab functions higher up the path can conflict with netlab

  % Set up vector of options for the Netlab EM optimiser.
  options = zeros(1,18);
  options(1) = 0;                 % This provides display of error values.
  options(5) = 1;                 % Reset singular covars to initial values
  options(14) = 100;              % Number of training cycles.

  Num_Centres = 5;

  mix = gmm(d, Num_Centres, 'full');
  mix = gmminit(mix, Data', options);
  mix = gmmem(mix, Data', options);

  Priors = mix.priors;
  Means = mix.centres';
  Covariances = mix.covars;  
end


switch lower(Covar_Type)
  case 'spherical'
    for i=1:Num_Centres
      eval(sprintf('Mean_%d = Means(:,%d);', i, i));
      eval(sprintf('Precision_%d = d/(sum(diag(Covariances(:,:,%d))));', i, i));
      eval(sprintf('Norm_%d = (Precision_%d/(2*pi))^(d/2);', i, i));
    end

  case 'diagonal'
    for i=1:Num_Centres
      eval(sprintf('Mean_%d = Means(:,%d);', i, i));
      eval(sprintf('Precision_%d = 1./(diag(Covariances(:,:,%d)));', i, i));
      eval(sprintf('Norm_%d = (prod(Precision_%d)/(2*pi))^(d/2);', i, i));
    end

  case 'full'
    for i=1:Num_Centres
      eval(sprintf('Mean_%d = Means(:,%d);', i, i));
      eval(sprintf('Precision_%d = inv(Covariances(:,:,%d));', i, i));
      eval(sprintf('Norm_%d = sqrt(det(Precision_%d))/(2*pi)^(d/2);', i, i));
    end

  otherwise
    disp(sprintf('Unknown covariance type %s', Covar_Type));
    exit(1);
end


